Vortices in Superfluid Fermi Gases Through the BEC to BCS Crossover 
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We have analyzed a single vortex at T — in a 3D superfluid atomic Fermi gas across a Feshbach 
resonance. On the BCS side, the order parameter varies on two scales: kp 1 and the coherence length 
f, while only variation on the scale of £ is seen away from the BCS limit. The circulating current 
has a peak value jmax which is a non-monotonic function of l/kpa s implying a maximum critical 
velocity ~ vf at unitarity. The number of fermionic bound states in the core decreases as we move 
from the BCS to BEC regime. Remarkably, a bound state branch persists even on the BEC side 
reflecting the composite nature of bosonic molecules. 



The recent discovery of vortices in 6 Li Fermi gases is 
a milestone in the development of fermion superfluid in 
atomic Fermi gases 0. For the first time, the phase co- 
herence of a superfluid atomic Fermi gas is demonstrated 
unambiguously. The atomic Fermi gases of alkali atoms, 
however, are very unusual many-body systems. Their in- 
teractions, which are related to s-wave scattering length, 
are highly tunable using Feshbach resonance. By tuning 
the inverse scattering length l/a s continuously from -co 
to +oo, the ground state of these systems changes from 
a weak coupling BCS superfluid to a molecular BEC 0- 
The region of infinite scattering (or a^ 1 =0) is particu- 
larly interesting, where the Fermi gas exhibits universal 
behavior, in the sense that the interaction energy scale 
at T = is independent of atomic details and is given 
by Fermi energy giving rise to a superfluid transi- 
tion temperature T c comparable to Fermi temperature 
Tp. In fact, it is found that at resonance, T c /ep ~ 0.20, 
the highest of all known fermion superfluids. 

Although all vortices of s-wave BCS superfluids have 
topologically invariant properties like quantized circula- 
tion of h/2M, other features like vortex core size, circu- 
lating current, and the bound state spectrum depend on 
dynamical details. It is natural to ask how the properties 
of a vortex change as one goes across the Feshbach reso- 
nance from the BCS to the BEC side, and how unitarity 
manifests itself in these properties. The purpose of this 
paper is to answer these questions using the Bogoliubov- 
deGennes (BdG) approach || in a three dimensional sys- 
tem for a wide Feshbach resonance 0|. We shall show 
that at T = 0: 

(A) On the BCS side, the order parameter A(p) for a sin- 
gle vortex exhibits two length scales: an initial rise on the 
scale of kp 1 Q for which we give an elementary analytical 
argument, and an eventual approach to its bulk value Ao 
on the coherence length scale £ = UvfJ Ao where vf is the 
Fermi velocity. At unitarity where Ao ~ £f = h 2 kp/2M, 
£ reduces to kp , and the two length scales coincide. On 
the BEC side the A(p) reaches its bulk value over the 
coherence length £ ~ 1/ ' -Jna s , where n is the density. 
We also find that the density n(p) is depleted in the vor- 
tex core, and its value n(0) at the center is dramatically 
reduced as one goes from the BCS to the BEC limit. 



(B) The circulating current j(p) around the vortex core 
has a peak value j max which is non-monotonic across the 
resonance and reaches a maximum precisely at unitar- 
ity. Its scale is set by the critical velocity v c which is 
determined by pair breaking on the BCS side, a single- 
particle effect, but by collective excitations on the BEC 
side. It is interesting to note that so far, j max is one of the 
very few properties of atomic Fermi gases we know that 
varies non-monotonically across resonance, in contrast to 
all thermodynamic properties which vary monotonically 
and smoothly. 

(C) We find that unitarity represents the most robust 
superfluid state in the entire BCS-BEC crossover. Not 
only does one obtain the highest T c at unitarity (which 
is, however, not too different from the T c value for 
all l/fepds > 0), but one also obtains the highest critical 
velocity v c ~ vf- 

(D) We find that as we go from the BCS to the BEC 
limit, the number of fermionic bound states in the vortex 
core decreases, with a corresponding increase in both the 
energy of the lowest bound state and their level spacing. 
Remarkably, we find that a bound state is observed way 
past unitarity, deep into the bosonic regime, which is 
unique to the molecular BEC. We also find that motion 
along the vortex core broadens a bound state of angular 
momentum i into a band with k z dispersion. 

The BdG approach was first used to study vortices in 
BCS superconductors in the classic work of Caroli, de- 
Gennes and Matricon see also Q. In the superfluid 
atomic Fermi gases, the BdG approach in two dimen- 
sions has been used in in the BCS limit and in [l2j |. 
The latter work is implicitly restricted to "narrow reso- 
nances" , which has very different physics from the wide 
resonances relevant to all current experiments 0, llflj . 
The vortex problem has also been studied using density 
functional |14| and hydrodynamic approaches, which 
are not microscopic and hence do not provide information 
such as (A) to (D) mentioned above. 

Bogoluibov-deGennes (BdG) approach: For a 

wide resonance, it is sufficient to use the single channel 
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model 0- The BdG equations in this case are 
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t A(r) \ / u n (v) 
A*(r) -f* J V u «( r ) 
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(1) 



where T — —h 2 V 2 /2M — p,, [i is the chemical potential, 
£" ra are the eigenvalues, u n and i> n are the eigenfunctions 
normalized by 



d 3 r [u* m (r)u n (r) + v* m (r)v n (r)} = 5 V 



(2) 



The order parameter A(r) and the chemical poten- 
tial p are determined by the self-consistency equation 
A(r) = gJ2 n Un ( r ) v n( r ) an d the average density n — 
2 Ji n I d 3 r\v n (r)\ 2 . The is restricted to < E n < 
E c where E c is a high energy cutoff (discussed below). 
^From now on we measure all energies in units of e F and 
lengths in units of kp 1 , so that the bare attraction g has 
units of kp/eF- The cutoff E c and the corresponding 
g(E c ) should be chosen to satisfy 
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so that the low energy effective interaction is described 
by the scattering length a s 10]. 

We work in cylindrical co-ordinates r = (p, z) in a 
gauge in which A(r) = A(p)e~ lB . Working in a cylindri- 
cal box of radius R and height L, our normalized wave- 
functions are of the form u n (r) — u n (p) e M e %kz z / V 2irL 
and v n (r) = v n (p)e t( > e+1 ^ e e tk: > z /V2irL so that decou- 
ples into different I and k z sectors We further ex- 
pand the radial wavefunctions u n (p) = Ylj c n j4>je(p) and 
v n (p) — Jlj d n j4>ji+i(p) m the orthonormal basis set 
(j)ji{p) = V2Je(aj£p/R)/[RJi + i(aji)] where a je is the 
j th zero of Ji(x). The BdG eigenvalue equation now re- 
duces to a matrix diagonalization problem 
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= E n 
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where TP 



,33 _ 



= {a 2 jR 2 + k 2 z - p)8jj> and A\ 
f dpA(p)(f)ji(p)(f)j/( + i(p) While the different i and k z 
sectors are completely decoupled in (@J, they are cou- 
pled through the self-consistency equations |l6(. Since 
the BdG equations are invariant under E n — > — E n , 
u n (r) — > v*(r), v n (r) — > — u*(r), one can get the pos- 
itive energy eigenfunctions for negative I by looking at 
the negative energy eigenfunctions for positive £. This 
reduces the computational effort by half. 

Ideally one would like to take E c — > oo and obtain so- 
lutions independent of this cutoff. In practice the size of 
Hilbert space grows like RL^/E- and to make the calcu- 
lation manageable we choose E c — 9e_F, R — 25k p 1 and 
L = 10k p 1 . We have checked that for this choice of cut- 
off, which is already larger than that in our results 



Q. 

< 



0.4 



1 .0 


------ ■ - " — ™ 




0.0 




"^-1.0 















k F p 



10 



15 



FIG. 1: The order parameter normalized by its value far from 
the vortex for three different couplings 1/kFa.s = —1,0,1. 
Note the clear emergence of two length scales in the BCS 
limit 1/kFds = — 1. 



in the uniform case are no more than 5% different from 
the infinite cutoff answers at unitarity. 

Vortex core structure: The order parameter A(p) 
is plotted in Fig. 1 for various values of l/kpa s . The 
weak oscillations in A(p), most prominent in the BCS 
regime, are likely finite size effects with no physical sig- 
nificance [TtI ]. Another very interesting aspect of the 
BCS regime, clearly visible in in fig. Q is the pres- 
ence of two length scales kp 1 and £ in the T = 
result for A(p). This is in marked contrast to the 
Ginzburg-Landau result where Aql(p) ~ tanh(p/£). 
While this effect was, in fact, recognized in the early 
superconductivity literature ||| in an Eilenberger calcu- 
lation, we provide an elementary analytical derivation 
here. Close to the origin, one may ignore A relative to 
the kinetic energy terms in the BdG equations and find 
U£ — Aj 1 Ji(k,Fp), vi — Aj 1 Ji + i{kpp) where Ai are con- 
stants. For small p we thus find A(p) = guo(p)v±(ft) = 
A^ 2 kFP- We next determine A§ as follows. For large 
p, A(p) = A and u ~ (l/^/B p) cos(k F p)exp{- p / '£) 
and vq = (1/\/Bqp) sm(kpp)exp(— p/£ ) (upto irrelevant 
phase shifts) where Bo is another constant. Match- 
ing the large and small p solutions at p = £, we get 
Aq ~ (Bokp 1 ) 1 / 2 0]. Finally using the normalization 
condition @ we can fix the constant Aq = £/(fcp£)^ 
which leads to the result A(p) ~ A^kpp for small p. 
Thus in the BCS limit, the initial slope of A(p) is set by 
kp 1 although the eventual approach to its uniform value 
A is on a second scale of £ = frvp/A ^> kp 1 . 

We see that outside of the BCS regime there appears 
to be a single length scale in A(p) as seen in Fig. 1. 
As the coupling increases toward unitarity the order pa- 
rameter value at large p, Ao increases toward cf and 
the scale £ shrinks to ~ kp 1 . The overall behavior of 
A(p) across the resonance is shown in Fig. EJa). At 
the same time the density profile n(p) around the vor- 
tex evolves as shown in Fig. Near the center of the 



3 



1.5 

Li. 

S 1 

< 

0.5 



(a) 



1.25 







"1/(k F a s ) = 1.0" 



0.5 



0.0 



0.5 



1.0 



k F p 




FIG. 2: (a) The order parameter profile and (b) the density 
profile for different values of l/{kpa s ). Inset: the density at 
the center of vortex as function of l/(fci?a s ). 
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FIG. 3: (a) The current distribution for three different val- 
ues of l/{kF<is) (b) The peak current jmax as function of 
l/(k F a a ). 



vortex n(p) ~ n(0) + ap 2 . The density at the center n(0) 
is a strongly decreasing function of l/kpa s as shown in 
the inset of fig. Gib) dropping from approximately 0.8n 
at l/kpa s = — 1, corresponding to a nearly "full" core, 
to O.ln at l/kpa s = +1, which is a nearly "empty" core. 

Circulating current: The current circulating around 
the vortex core is j = p s v s , where p s ~ A 2 , and v s = 
(h/2Mp)9. Far away from the vortex core, j ~ 1/p since 
A — > Aq. On the other hand, near the center of the 
vortex, A ~ p, so that j ~ p. This implies that the 
current must have a peak (j max ) at a distance p* from 
the center, and p* can be taken as a measure of the size of 
the vortex core. Intuitively, as the current increases with 
decreasing p it eventually reaches the critical current at 
p* , so that for p < p* the superfluid order parameter is 
destroyed. 

Formally the current density is given by 



rdk z j 



(5) 



The profiles of for various l/kpa s are shown in 
Fig. EI a). They all show a peak, as explained above. We 
find that the location p* of the maximum current shows a 
weak non-monotonic behavior in the vicinity of unitarity, 
which is roughly consistent with the non-monotonic be- 
havior of the coherence length £ as a function of 1/ fcj?a s 
predicted earlier; see Fig. 3 of [ijj. Details will be pub- 
lished separately |2cj . 

The peak current j rn ax, however, shows a non- 
monotonic behavior as a function of l/kpa 3 , peaking at 
unitarity as shown in Fig.[3{b). We may understand this 
non-monotonic behavior as follows. On the BCS side, 
the critical velocity, determined by depairing, is A /kp 
and thus j ma x ~ nvp(Ao/ep). The critical current thus 
increases as one moves toward unitarity from the BCS 
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FIG. 4: The BdG spectrum in the k z = sector as a function 
of I for (a) l/(fejro s ) = 0.0, i.e., unitarity and (b) l/(fcFO s ) = 
1.0, i.e., on BEC side, (c) The minigap [t — k z = 0) as a 
function of l/(kpa s ). (d) The k z dependence of the bound 
state energy for I = at unitarity. 



side as Ao increases. At unitarity, universality dictates 
that the critical velocity must be of order vp so that 
jmax ~ nvp. As one goes toward the BEC side the crit- 
ical current is now determined not by pair breaking but 
rather by the collective excitations. The Landau crite- 
rion suggests v c ~ vp\/kpa s which leads to a j max which 
decreases with increasing l/kpa s . We thus find the in- 
teresting result that the mechanism for destruction of 
superfluidity, as reflected in the maximum current j ma x > 
changes as one goes across the resonance from pair break- 
ing on the BCS side to collective on the BEC side. 

Spectrum: We now turn to the fermionic bound 
states in the vortex core and their evolution through the 
BCS-BEC crossover. In the BCS limit the results are 
very well known from CDM They showed that for 
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each £ and k z there is a bound state, with energy less 
than Ao, in the core of the vortex. The lowest energy 
fcrmionic excitation (with quantum numbers £ — and 
k z = 0) has a "minigap" Ag/2e F <C A . As we increase 
l/kpds we see some changes in the spectrum. Our re- 
sults at unitarity l/kpCbg = and in the BEC regime 
l/fc_pa s = 1 are shown in Fig. 01(a) and (b) respectively. 
For clarity we show the spectrum only as a function of 
angular momentum £ and at fixed k z = (fc 2 -dispersion 
is discussed later). 

As long as the chemical potential p > (which includes 
unitarity), the bound (continuum) states are those with 
energies smaller (larger) than A . The bound state spec- 
trum at unitarity (Fig.0Ja)) is not qualitatively different 
from the BCS limit, except that both the minigap and 
the level spacing are larger, and therefore one has fewer 
bound states. Remarkably, the minigap continues to fol- 
low Aq/2cf even through unitarity as shown in Fig.0Jc). 

Once the chemical potential p < 0, as in the BEC 
regime shown in Fig. EJb), the continuum of fcrmionic 
excitations exists for E > (\p\ 2 + A 2 ,) 2 |2j. We can still 
define a fermionic bound state by demanding that the 
corresponding wavefunction decays exponentially to zero 
away from the vortex core. If such states exist, then it 
can be easily shown that their energy must lie in the in- 
terval \p\ < E < (\p\ 2 + A^)? . Remarkably, we find such 
a bound state well into the BEC regime as evident from 
Fig. EJb). It is amusing that the off-diagonal potential 
A(p) produces an Andreev bound state even at an energy 
larger than the maximum value of the potential Ao, but 
we must emphasize that this bound state is below the 
continuum which starts at (\p\ 2 + A 2 ,) 5 . The fermionic 
bound states are a unique consequence of composite na- 
ture of bosons in the molecular BEC and are absent in 
atomic BEC. 

As a result of motion along the vortex axis, each 
state for fixed angular momentum £ shown in Fig. 0] (a) 
and (b) actually broadens into an energy band with k z - 
dispersion. Our calculated k z dependence of the bound 
state energy for the I = bound state at unitarity is 
shown in Fig.QJd). The energies (which are discrete due 
to the finite L in the z-direction) continue to follow the 
Caroli et. al. BCS-limit prediction E /(l - k 2 )i even 
at unitarity . 

Finally, we note that there is a branch of bound states 
(shown with smaller dots in Figs.QJa) Bib)) which lies 
below the fermionic energy gap in the bulk but differs 
from the vortex core states (which exist only for £ > 0) 
in that it exists both for positive and negative £. From 
their wavefunctions we deduce that these states are not 
related to the vortex core, but are in fact trapped near 
p = R due to the suppression of the order parameter 
by the hard wall boundary condition ji^]. The details 
of these as well as the vortex bound state wavefunctions 
will be described separately [20(. 



Conclusions : The physics of a vortices in the BEC to 
BCS crossover has led to interesting results (A) through 
(D) summarized in the introduction. The study of the 
order parameter profile in the vortex core and the cir- 
culating current may require new experimental methods 
which will allow more detailed diagnostics. The existence 
of fermionic bound states should be detectable by spec- 
troscopic means, or through the damping of sound, or 
through effects of dissipation in vortex dynamics. We 
hope that the unique properties of vortices that we point 
out here will stimulate further experimental and theoret- 
ical studies on new methods to probe strongly interacting 
Fermi gases across Feshbach resonance. 

This work is supported by NASA GRANT-N AG8- 1 765 
and NSF Grant DMR-0426149. 
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